Combined effects of wind, air temperature and snowfall on snow accumulation in a subalpine forest
Introduction
Forest cover extends over approximately half of North America’s snow-covered region, governing snowpack accumulation and ablation, thereby influencing the timing and magnitude of runoff generation from snowmelt.
There is a need for robust models of snow redistribution by vegetation and wind to estimate snow accumulation in mountain forests. To achieve this, a comprehensive understanding of snow redistribution processes is required. However, existing snow interception parameterizations are based on limited observations in distinct climates and forest structures. Rapid changes in climate and forest ecology illustrate the pressing need to assess whether existing snow interception parameterizations are suitable for predicting snow accumulation in diverse and changing environments.
Intercepted snow in the canopy is subjected to higher rates of sublimation compared to subcanopy snow due to greater surface area, higher wind speed, and solar exposure (Pomeroy et al., 1998). Across the Northern Hemisphere, researchers estimate that 25 to 45% of annual snowfall may be lost to the sublimation of intercepted snow from the canopy (Essery et al., 2003). Correctly determining the fraction of snowfall intercepted in the canopy is crucial for estimating interception losses by sublimation (Pomeroy et al., 1998). In addition, the time that snow resides in the canopy and is subject to sublimation is dependent on rates of unloading, melt, drip, and resuspension of snow (Hedstrom & Pomeroy, 1998; Katsushima et al., 2023; Lumbrazo et al., 2022; Storck et al., 2002).
The theory underpinning current snow interception parameterizations is based on observations ranging from warm maritime (Andreadis et al., 2009; Storck et al., 2002) and cold continental (Ellis et al., 2010; Hedstrom & Pomeroy, 1998; Roesch et al., 2001; Satterlund & Haupt, 1967) climates generally characterized by dense forest canopy. Accurate simulations of forest snow accumulation have been achieved if the parameterizations are applied in similar climates to where they were developed (Lundquist et al., 2021; Rasouli et al., 2019; Roth & Nolin, 2019) or if they are combined into a hybrid parameterization and assessed at global and regional scales in a wide range of climates (Essery et al., 2003; Gelfan et al., 2004). Although accurate performance has been achieved across different climates in some studies (Essery & Pomeroy, 2004; Gelfan et al., 2004), other snow model comparisons (Krinner et al., 2018; Rutter et al., 2009) have shown reduced performance. The decision in earth system models to use parameterizations derived in warm or cold climates is often based on a simple temperature-based step function (Essery et al., 2003; Gelfan et al., 2004) and may require modification to better represent more transitional climates and forest types. The original theory has also been simplified over time, i.e. the increase in canopy coverage with increasing wind speed is not included in more recent parameterizations Roth & Nolin (2019). Updates by Gelfan et al. (2004) to combine the Hedstrom & Pomeroy (1998) and Storck et al. (2002) parameterizations is not typically utilized in recent studies (Krinner et al., 2018; Rutter et al., 2009). The omission or simplified representation of processes and reliance on empirical calibrations likely contribute to model uncertainty when applied in climates and forests where other processes become important (Krinner et al., 2018; Lumbrazo et al., 2022; Lundquist et al., 2021; Moeser et al., 2015; Roth & Nolin, 2019; Rutter et al., 2009).
Figure 1 shows the difference in the change in interception efficiency across different snowfall event sizes for three common models. Interception efficiency (interception/snowfall) declines with increasing snow load (Hedstrom & Pomeroy, 1998; Storck et al., 2002) or initially increases and then is followed by a decline (Moeser et al., 2015). The underlying theory of the Moeser et al. (2015) parameterization stems from the Satterlund & Haupt (1967) study who observed an initial increase in the rate of intercepted snow, as snowflakes bridge gaps between needles. It may also be inferred that during the small near 0°C snowfall events observed in Satterlund & Haupt (1967) the majority of snow may have melted immediately due to a warm canopy resulting in low interception efficiency. The initial low interception efficiency was followed by an increase and then flattening off of the interception rate as branches bend due to the weight of snow which Satterlund & Haupt (1967) represented by a numerical analytical sigmoidal function. In the observations by Hedstrom & Pomeroy (1998), snow interception efficiency starts high and then declines. Hedstrom & Pomeroy (1998) hypothesis the shape of this curve is due to a decrease in canopy contact area and change in the incoming snowfall angle of impact as branches bend downward. This relationship was best represented using an inverse exponential function shown in Figure 1. The Hedstrom & Pomeroy (1998) parameterization therefore differs from the (Moeser et al., 2015; Satterlund & Haupt, 1967) sigmoidal function as it does not include a representation for the initially slow interception rate. Andreadis et al. (2009) developed a snow interception model using data collected by Storck et al. (2002) in dense old growth forest in the maritime climate of southwestern Oregon, USA. This method builds off the maximum canopy snow load theory proposed in Hedstrom & Pomeroy (1998) but makes additional modifications to include a step function based on temperature. Here, snow interception efficiency, was found equal to a constant of 0.6 based on snow interception observations from Storck et al. (2002) in southern Oregon.
Maximum interception capacity decreases (Hedstrom & Pomeroy, 1998) or increases (Storck et al., 2002) with increasing air temperature (Figure 2). Hedstrom & Pomeroy (1998) proposed that fresh snow density, which may be described as a function of air temperature (Hedstrom & Pomeroy, 1998), plays an important role in governing the interception capacity. Storck et al. (2002) limit, \(L\) as being less than or equal to the maximum interception storage \(L_{max}\) using a step function of temperature based on observations of warmer snow having more cohesion to the canopy.
More recent work by Katsushima et al. (2023) collected measurements of snow interception using a weighed tree for a warm-humid coastal environment in Japan. They observed a decline in interception efficiency with increasing wind speed, they attributed to increased hydrometeor velocity and bouncing on impact. While not mentioned in this study, the decrease in interception efficiency may also be due to wind induced unloading. They did not observed a maximum interception capacity within their measurement range of 0-25 mm. Although for temperatures above 0 they could see a decline in interception efficiency above 10 mm maybe due to branch bending + melt rates. Katsushima et al. (2023) suggest air temperature and wind speed alone are insufficient to describe interception efficiency and hypothesize that particle shape may be an improved predictor but did not have the observations to test this. Some of the limited model performance reported by Katsushima et al. (2023) may be attributed to a result of their interception measurements including unloading due to melt and wind.
Previous studies have collected measurements of interception efficiency over snowfall events ranging from hourly (Storck et al., 2002) to weekly timesteps (Hedstrom & Pomeroy, 1998). The different measurement time intervals vary in the amount of time possible for ablative processes and may influence model estimates of interception. As a result, some of the interception measurements inevitably include some amount of ablation. Despite the inclusion of unloading in the interception parameterizations developed in these studies, they are often combined with additional unloading parameterization in earth system models (Clark et al., 2015; Ellis et al., 2013) leading to some potential of double counting of the ablation process.
Uncertainties also arise in the scaling of point or branch scale measurements to the plot scale (Staines & Pomeroy, 2023).
Forest structure governs the interception efficiency observed at a given location (Hedstrom & Pomeroy, 1998; Roth & Nolin, 2019). Metrics used in common snow interception parameterizations (Hedstrom & Pomeroy, 1998; Storck et al., 2002) to describe forest structure include canopy cover and leaf area index (LAI). Leaf area index is defined by Chen et al. (1997) as one half the total green leaf area per unit ground surface area. Canopy cover is defined in Hedstrom & Pomeroy (1998) as the fraction of sky not visible by the instrument from under the canopy. While more detailed forest structure metrics exist derived from detailed LiDAR scans (Helbig et al., 2020; Roth & Nolin, 2019), often they are not available at regional extents required to run hydrological models.
Add sentence on Staines & Pomeroy (2023).
Several processes govern the accumulation of snow in mountain forests, and the importance of individual processes may differ depending on climate and forest structure (Gelfan et al., 2004; Hedstrom & Pomeroy, 1998; Moeser et al., 2015; Staines & Pomeroy, 2023). Therefore determining the dominant processes in varying climate and forests structures is important to guide model decision makers on what existing model parameterizations to choose or if a hybrid approach may be appropriate.
New observations of snow interception and ablation processes will help determine if existing theories are applicable in differing climates and diverse forest structures.
The novelty of this study is the study site location in a windswept discontinuous subalpine ridge forest and the use of high temporal frequency automated measurements and discrete high spatial resolution measurements using aerial LiDAR to attempt to separate out interception from ablative processes. Timelapse cameras were also used to confirm absence of unloading during interception periods.
Exisiting theory relies on the green smear approach or using LAI + CC which does not appropriately describe the heterogeneity of discontinuous canopies
Objective: To assess the influence of canopy structure and meteorology on snow interception processes in a windswept subalpine forest.
Research Questions:
- Are the theories and assumptions of existing snow interception parameterizations true for field measurements collected in a continental subalpine forest?
- What are the dominant processes that control snow interception in a subalpine forest?
Specific Questions
Is meteorology important for governing I/P?
How does hydrometeor trajectory angle influence apparent canopy structure?
How does the association between canopy structure and I/P vary across the hemisphere for an event with moderate wind speed?
is this association stronger when there is snow in the canopy i.e., branches have been bridged or reduced due to branch bending?
Theory
Figure 3 shows the mass exchange of snowfall between the forest canopy and the snowpack on the ground surface. The storage of snow water equivalent (SWE) is represented here with respect to the canopy (\(L\), mm) or the surface snowpack (\(S\), mm). Fluxes that are repeated between the canopy and snowpack control volume have a superscript to specify what control volume they refer to (i.e. \(q^{veg}\) refers to the vegetation control volume and \(q^{snow}\) refers to the surface snowpack control volume).
The change in canopy SWE storage, \(L\) (mm), may be represented as:
\[ \frac{dL}{dt} = q_{sf} - q_{tf}(L) - q_{unld}(L) - q_{drip}(L) - q_{wind}^{veg}(L) - q_{sub}^{veg}(L) \tag{1}\]
where \(q_{sf}\) (mm s-1) is the above canopy snowfall rate, \(q_{tf}\) (mm s-1) is the throughfall rate, \(q_{wind}^{veg}\) (mm s-1) is the wind transport rate in our out of the control volume (typically assumed to be negligible in the literature), \(q_{sub}^{veg}\) (mm s-1) is the intercepted snow sublimation rate, \(q_{unld}\) (mm s-1) is the canopy snow unloading rate and \(q_{drip}\) (mm s-1) is the canopy snow drip rate due to canopy snowmelt. \(q_{wind}^{veg}\) and \(q_{sub}^{veg}\) may be a positive or negative flux. Where all of the above rates are a function of snow load (\(L\)), which is how much snow is present in the canopy.
The trajectory angle, \(\theta_h\) of a hydrometeor as the departure in degrees (°) from a horizontal plane (i.e., -90° for vertical snowfall), may be calculated as:
\[ \theta_h = \arctan \left(\frac{-v_h(D_h)}{x_h(u_z)}\right)*\frac{180}{\pi} \tag{2}\]
where \(v_h\) is the terminal fall velocity of the hydrometeor (m s-1), which is a function of the hydrometeor diameter, \(D_h\), \(x_h\) is the horizontal change in the hydrometeor (m s-1) which is a function within canopy wind speed, \(u_z\) at height above ground, \(z\).
Above the top of the canopy, wind flow is approximately logarithmic with height. Cionco (1965) show that, \(u_z\) may be approximated using the exponential formula:
\[ u_z = u(u^*,z,z_o,d)\cdot exp\left[a\cdot\left(\frac{z_c}{h_c}-1\right)\right] \tag{3}\] where \(u\) is the horizontal wind speed at the top of the canopy which is a function of the friction velocity \(u^*\), height above ground, \(z\), roughness length, \(z_o\), and the displacement height of the canopy, \(d\). \(a\) is an attenuation coefficient that increases with increasing leaf area and decreases as the mean distance between individual trees increases, \(z_c\) is the height above ground of \(u_z\), and \(h_c\) is the average height of the canopy elements. Cionco (1972) suggest values for \(a\) of 1.01 for small needleleaf trees, Zhu et al. (2001) provide methods to calculated \(a\) based on canopy density and Parviainen & Pomeroy (2000) provide a method to calculate \(a\) using observations from two jack pine (Pinus banksiana) stands. Figure 4 shows the increase in trajectory angle calculated using Equation 2 with a constant hydrometeor velocity of 1 m s-1 and a horizontal velocity equal to mid-canopy wind speeds ranging from 0-20 m s-1.
Methods
Study Site
This study was conducted at forest plots surrounding the Powerline Station (PWL) and Forest Tower Station (FT) located within Fortress Mountain Research Basin (FMRB), AB ~ 51 ◦ N, 2100 m asl. (Figure 5). The species of tree surrounding the PWL and FT stations include sparsely spaced coexisting subalpine fir (Abies lasiocarpa) and Engelmann spruce (Picea engelmannii), with a proportion of 70% and 30% respectively (Langs et al., 2020). In the early 1900s the majority of the FMRB vegetation burned during a large forest fire that affected most of the Kananaskis Valley. Following the fire, the forest surrounding the PWL and FT stations has naturally regenerated with minimal disturbance.
Meterological Measurements
A flux tower at the FT Station provided measurements of air temperature (Campbell Scientific HMP155A), relative humidity (Campbell Scientific HMP155A), wind speed (Campbell Scientific Ultrasonic-86000 & Metone 014A) and direction (Campbell Scientific Ultrasonic-86000) at 15-min time intervals 4.3 m above the ground (Figure 5 (b)). The PWL station, located 120 m to the northwest in a forest clearing shown in Figure 5 (b) provided measurements of snowfall rate (\(q_{sf}\)), additional wind speed measurements and hydrometeor size and velocity. Snowfall was measured using a weighing precipitation gauge 2.6 m above ground (OTT Pluvio) corrected for undercatch following phase correction by Harder & Pomeroy (2013) and catch efficiency by Macdonald & Pomeroy (2007). A Metone 014A anemometer provided wind speed measurements at a height of 2.6 m, while the RM Young anemometer recorded wind speed and wind direction at a height of 5.2 m. A disdrometor (OTT Parsivel2) provided measurements of hydrometeor particle size and velocity.
Lysimeter Data
Three subcanopy lysimeters (SCL) installed surrounding the FT Station provided fifteen minute interval measurements of throughfall plus unloading (SCL locations are shown in Figure 5 (b)). For select events where ablative processes, \(q_{unld}(L)\), \(q_{drip}\) and \(q_{wind}^{veg}\) could be considered negligible, the subcanopy lysimeters were inferred to provide measurements of throughfall. Timelapse imagery, a weighed tree and in-situ observations were used to ensure ablation of snow intercepted in the canopy or snow on the ground was minimal over each of the selected events. The SCLs consisted of a plastic bucket with an opening of 0.9 m 2 and depth of 20 cm (e.g., Figure 6 (a)) suspended from a load cell (Intertechnology 9363-D3-75-20T1) attached to an aluminum pipe connected between two trees. The load cell which measures kilograms was scaled to kg m-2 by dividing by the cross-sectional area of the SCL opening. The canopy structure surrounding three SCLs and was measured using hemispherical photography (Nikon Coolpix 4500 and EC-F8 hemispherical lens) and the hemispheR R package Chianucci & Macek (2023) and is shown in Table 1.
A weighed tree lysimeter, shown in Figure 6 (b) and Figure 6 (c), measured the weight of canopy snow load, \(L_{wt}\) (kg). A live subalpine fir (Abies lasiocarpa) tree was cut and suspended from a load cell (Artech S-Type 20210-100) at the beginning of the 2022 and 2023 water years which recorded the weight of the tree. The bottom of the tree was sealed to limit some transpiration and to prevent spinning and abrupt impacts of the free-hanging tree, the base of the tree was attached to a support system that allows for vertical movement but limits abrupt horizontal movements. The weight of snow in the weighed tree in kg was scaled to an areal estimate of canopy snow load (\(L\), kg m-2) using Equation 1 and measurements of areal throughfall (kg m-2) from manual snow surveys and snowfall from the PWL Station snowfall gauge (see description of method in Hedstrom & Pomeroy, 1998).
UAV Data Processing
Two uncrewed aerial vehicle (UAV) lidar surveys were conducted before and after the snowfall event to be used for the measurement of snow accumulation and computation of canopy structure metrics. The UAV (FreeFly Alta X) was equipped with a REIGL miniVUX-2 airborne laser scanner payload, an Applanix APX-20 inertial measurement unit (IMU) and global navigation satellite system (GNSS). The UAV was flown 90 m above the ground at a speed of 3 m s-1 following a preprogrammed flight trajectory shown in Figure 5 (b).
The REIGL miniVUX-2 laser operates at a near infrared wavelength with a laser beam footprint of 0.160 m x 0.05 mm (at 100 m above ground). The accuracy and precision of the miniVUX-2 is described by REIGL for a lab environment of 0.015 m and 0.01 m respectively (at 50 m above ground). The miniVUX-2 was configured with a laser pulse repetition rate of 200 kHz, field of view of 360°, scan speed of 31.09 revolutions s-1 and an angular step width of 0.0558°, resulting in an expected an average point cloud density of 107 returns m-2 for each flight path.
Georeferenced point clouds with x, y, and z coordinates for each laser return were generated following methods outlined by Harder et al. (2020) and Staines & Pomeroy (2023) to reconcile survey lidar, IMU and GNSS data. A ground based GNSS system was positioned on a permanent monument during each survey and underwent precise point positioning (PPP) correction by Natural Resources Canada (2024). Differential GNSS correction of the UAV trajectory was conducted using the ground based PPP GNSS observations and the POSPac UAV software. The UAV-lidar point clouds were then transformed from a sensor referenced coordinate system to a georeferenced coordinate system (EPSG:32611 - WGS 84 / UTM zone 11N) using the RIEGL Riprocess Software. A consistent vertical offset of up to 6 cm between UAV-lidar flight lines was observed in the resulting point clouds on March 13th and 14th, 2024 and was attributed to IMU drift. This offset between flight lines was corrected using the BayesMap stripalign software, which adjusts the vertical elevation of the point cloud to the ground control points (GCP) collected across the study site using a differential GNSS rover. After strip alignment, the mean elevation bias (lidar minus GCP) was 0.000 m and the RMS error changed from 0.055 m to 0.038 m March 13th and changed from 0.033 m to 0.029 m on March 14th. The point cloud density ranged from ~1200 returns m2 in open clearings to ~2200 m2 in sparse forest for both the March 13 and 14th surveys after all flight paths were combined.
Quality control, ground classification and calculation of the change in between two UAV-lidar point clouds was conducted using the LAStools software package (LAStools, 2024). Point clouds were clipped to a polygon with a 20 m buffer bounding the UAV transect shown in Figure 5 (b). The ground classification was conducted using the “lasground_new” function (LAStools, 2024) for both the pre and post snowfall event point clouds, with a step size set to 2 m and 8 substeps (ultra_fine setting). The offset and spike options were set to remove points that are more than 0.1 m above or below the initial ground surface estimate surface which “lasground_new” fits to the last returns. This function is based on an algorithm outlined by Axelsson (2000), describing the process of making the initial ground surface element.
Snow Surveys
In-situ Snow Depth and Density
Manual snow surveys provided measurements of subcanopy throughfall over the forest surrounding FT Station after a snowfall event. These throughfall measurements were used to scale the weighed tree and for accuracy assessment of the UAV-lidar fresh snow depth measurements. If a pre event crust layer was present the depth of post event fresh snow accumulation above the crust layer were interpreted as throughfall over the event. In the absence of a defined crust layer, the difference in pre and post event snow depth to ground was interpreted as event throughfall. A 1000 cm3 snow density wedge sampler (RIP Cutter, https://snowmetrics.com/shop/rip-1-cutter-1000-cc/) was used to measure the density of the fresh snow layer, \(\overline{\rho_{tf}}\) (kg m-3). The throughfall depth measurements, \(\Delta HS\) were converted to snow water equivalent (SWE) using the following equation:
\[ \Delta SWE_{tf} = \Delta HS \cdot \overline{\rho_{tf}} \tag{4}\]
Fourteen snow surveys (seven pre and post snowfall event pairs) that had minimal evidence of ablation were select for calibration of the weighed tree. An average of 40 snow depth samples were taken for each snow survey following two 60 m linear transects north and south of the FT Station. The forest canopy surrounding the snow survey stations had LAI ranging from 0.21 – 3.91 (-) and canopy coverage ranging from 0.14 – 0.86 (-) measured using hemispherical photography.
UAV-Lidar Snow Depth
The change in elevation between the two UAV-lidar surveys was interpreted as the increase in snow accumulation, \(\Delta HS\) over the snowfall event. This change was calculated using a point-to-grid subtraction method, similar to the approach used by Deems et al. (2013) and Staines & Pomeroy (2023), using the “lasheight” function from the LAStools (2024) software. The pre snowfall event point cloud from “lasground_new” by “lasheight” to construct a “ground” TIN. Subsequently, the height of each post snowfall event point above the ground TIN, resulting in a point cloud representing \(\Delta HS\). This point cloud was then converted into a raster of \(\Delta HS\) with a grid cell resolution of 5 x 5 cm using the “las2dem” function. Further quality control and resampling of the 5 cm raster of \(\Delta HS\) was conducted using the R package Terra Hijmans (2024). Regions that were disturbed over the snowfall event during the in-situ snow survey and values that exceeded the .999th quantile were removed. To help remove any remaining noise a 25 cm \(\Delta HS\) raster was generated by computing the median of the 5 cm \(\Delta HS\) values within each 25 cm grid cell.
Snow Interception
Throughfall measurements collected form manual snow surveys (\(q_{tf}\cdot \Delta t\)) and the subcanopy lysimeters (\(q_{tf}\)) were used to estimate the amount of snow intercepted in the canopy. During calm snowfall periods Equation 1 was be simplified to estimate the amount of snow intercepted in the canopy:
\[ \frac{dL}{dt} = q_{sf}-q_{tf}(L) \tag{5}\]
This method was preferred, compared to measurements from the weighed tree lysimeter, as the subcanopy lysimeters were not influenced by sublimation losses from snow intercepted in the canopy.
Interception efficiency, \(\frac{I}{P}\) (-), which is the fraction of snow intercepted over a discrete time interval, \(\Delta t\) was calculated as:
\[ \frac{I}{P} = \frac{\Delta L}{\overline{q_{sf}} \Delta t} \tag{6}\]
where \(\Delta L\) (mm) is the increase in canopy load over a discrete time interval and \(\overline{q_{sf}}\) (mm), is the average snowfall rate over \(\Delta t\).
Results
The influence of meteorology on snow interception
The meteorology, state of canopy snow load, and interception efficiency was measured for 26 distinct snowfall events. The histograms presented in Figure 7, portray the frequency distribution of 15-minute measurements of air temperature, relative humidity and wind speed from the FT Station, snowfall rate from PWL Station, and interception measurements calculated using SCL 1-3 and snowfall gauge at PWL Station during all periods of snowfall during the 26 snowfall events (See station locations in Figure 5 (b)). The average interception efficiency observed across all periods of snowfall was 0.64 and ranged from 0 to 1. The air temperatures spanned from -25°C to 1°C, with a predominant concentration of observations within the -10°C to 0°C range. Corresponding wind speed values ranged from 0 to 4.5 m s-1 with a mean of 1.2 m s-1.
The state of canopy storage with event cumulative snowfall, measured using SCL 1-3 and snowfall gauge, is shown in Figure 8 (a) for each of the 26 snowfall events. Canopy storage was observed to increase linearly with increasing snowfall without evidence of reaching a maximum canopy storage capacity. Variation in the slope of each line in Figure 8 (a), is attributed to differences in the canopy structure surrounding the SCL instruments and differences in meteorology and antecedent canopy snow load between the individual events. The SCL with higher canopy coverage was expected to have consistently higher canopy storage values. However, as depicted in Figure 8 (a) the largest event cumulative snowfall of ~ 50 mm, shows the lysimeter with the highest canopy coverage (0.82) exhibited a canopy storage value slightly less than the lysimeter with canopy coverage of 0.78. This difference between the subcanopy lysimeters, which is inconsistent with their individual forest structure, is likely result of event wind speed and the direction which could alter the apparent forest structure of each lysimeter location. In addition, the initial canopy snow load may have been different between the locations. The absence of canopy storage and interception efficiency measurements for certain troughs during specific events was caused by damage to the subcanopy lysimeter wiring due to animal interference and snow load.
For each of the 26 snowfall events the average interception efficiency was calculated and is shown with the corresponding event snowfall total in Figure 8 (b). The average interception efficiency for SCL 1-3, across all 26 events, is shown in Figure 8 (b) to increase proportionally with the canopy coverage of each lysimeter. However, looking at the individual events, interception efficiency observed at each lysimeter does not always align with the respective canopy coverage value as was also shown in Figure 8 (a). This suggests additional processes other than canopy structure are influencing subcanopy snow accumulation. For events with snowfall totals ranging from 0 to 5 mm, the average interception efficiency was 0.58. This increased to an average interception efficiency of 0.60 for events with snowfall totals between 5 and 10 mm, then decreased back to 0.58 for events with snowfall totals between 20 and 25 mm.
The average interception efficiency was also calculated over each 15-minute interval within each of the 26 snowfall events and was utilized to assess the association with meteorological variables, the state of canopy snow load and hydrometeor characteristics for the corresponding time interval. Air temperature and relative humidity observed at FT Station was observed to have little influence over interception efficiency measured using the SCL and snowfall gauge (Figure 9, A & B). Mid canopy wind speed at FT Station was observed to increase the mean interception efficiency slightly from 0.56 to 0.62 between wind speed bins of 0.25 and 1.25 m s-1 (Figure 9, C). This is thought to be due to an associated increase in canopy contact area as hydrometeor trajectory becomes more horizontal with increasing wind speed. The mean interception efficiency decreases wind speeds above 2 m s-1 to a minimum of 0.48 for the 3.75 m s-1 wind speed bin. The initial canopy snow load had a positive association with interception efficiency as the canopy filled with snow, increasing the mean interception efficiency from 0.57 for snow loads near 0 mm to 0.66 for snow loads around 7 mm. This was followed by a gradual decline in interception efficiency for snow loads greater than 7 mm to a minimum of 0.48 at snow loads above 18 mm (Figure 9, C). Hydrometeor diameter and hydrometeor velocity measured using the disdrometer at PWL Station did not have a strong association with interception efficiency (Figure 9, E & F).
The influence of forest structure on snow interception
Snow accumulation was measured over the UAV-lidar flight extent (Figure 5 (b)) after a 24 hour snowfall event with steady southerly winds. The meteorology from this snowfall event, totalling 28.4 mm between lidar flights on March 13th and 14th, is shown in Figure 10. The event was characterised by a transition from low rates of snowfall and air temperature near 0°C to higher rates of snowfall late afternoon on March 13 and coincided with air temperatures around -2.5 °C. Wind speed observed while it was snowing (\(q_{sf}\) > 0.1 m s-1) was predominately lower at FT Station compared to PWL Station with event averages of 1.27 and 1.75 respectively. Wind direction observed during snowfall was consistently from the south at both FT Station and PWL Station with an average of 193°. This snowfall event resulted in a range of throughfall depths measured across the study site, shown by the frequency distribution of throughfall in Figure 11.
The spatial distribution in I/P across the study site, calculated using Equation 6 and \(\Delta L\) from UAV-lidar measurements and \(\Delta sf\) from the Pluvio snowfall gauge is shown in Figure 12. Greater I/P is observed on the north (lee) side of individual trees which is interpreted to be a result of non vertical hydrometeor trajectories from the steady southerly winds observed over this event. This effect is more apparent within the PWL forest plot (Figure 12 (c)), compared to the FT forest plot (Figure 12 (d)). While average wind speed over the event was observed to be higher at the PWL Station compared to the FT Station, the position of the anemometer in a forest clearing and slightly higher off the ground likely contributed to the higher wind speed. Wind speeds within the PWL forest plot canopy are expected to be lower compared to the FT forest plot due to the higher canopy coverage and taller canopy height. The height of the canopy measured using UAV-lidar, is shown in Figure 12 (a) and fig-lidar-ch-ft for the PWL and FT forest plots respectively. The average height of the canopy surrounding the plot to the east of the PWL station is 10.51 and surrounding the forest plot around the FT Station is 7.12. Based on the large difference in forest structure it is interpreted that the differences in I/P between the two plots are primarily a result of differences in forest structure rather than localised wind flow patterns.
To determine how forest structure was associated with interception efficiency at different azimuth and zenith angles over the March 13-14 snowfall event, each portion of the hemisphere at each grid location was considered. The relationship between I/P and canopy contact number was found to be non-linear and thus the Spearman’s Correlation Coefficient, \(\rho_s\) was calculated to quantify the association between a single raster of I/P and 32’760 rasters containing the canopy contact number hemisphere for each portion of the hemisphere (azimuth [0, 1, …, 359], zenith angle [0, 1, …, 90]). A hemisphere plot was generated for both the FT and PWL forest plots using the \(\rho_s\) value from each azimuth and zenith angle combination and is shown in Figure 13. The highest correlation between canopy contact number and interception efficiency was found to the south between zenith angles of approximately 10 to 30° (trajectory angle -80 to -60°) for both PWL and FT plots. The high \(\rho_s\) around similar zenith angles between PWL and FT also provide evidence that wind speeds were similar between the two plots resulting in similar hydrometeor trajectory angles. The predicted hydrometeor trajectory angle of the event, calculated using Equation 2, was -35.32° using measured wind speed (4.5 m above ground) at FT Station and hydrometeor velocity of 0.9 m s-1 from PWL Station. The predicted hydrometeor trajectory angle is 25° to 45° more horizontal compared to the area of the hemisphere found to have the highest correlation between canopy contact number and I/P. This may be partially explained by three factors, 1) extrapolating wind speed measured within sparse canopy (FT Station) to regions characterised by higher stem density (PWL Station), where wind speed was likely lower leading to more vertical trajectory angles, 2) The conical shape of needleleaf trees within the PWL and FT forest plots, resulting in a larger surface area closer to the ground, and thus intercepted more snow closer to the tree an thus leading to higher correlation between I/P and contact number from zenith angles closer to vertical, 3) the theory proposed in Equation 2 assumes a linear trajectory, and does not consider non-linear patterns such as wind flow wrapping around tree elements, turbulent flow, or differences in wind speed with height.
Combined Effects of Meteorology and Forest Structure
The mean number of canopy contacts calculated for each forest plot for different hydrometeor trajectory angles [-90, -85, …, 0] averaged across each azimuth angle [0, 1, …, 359] and across all plot grid cells is shown in Figure 14 (a) to increase exponentially with increasing horizontal hydrometeor trajectory angle (estimated as a function of wind speed using Equation 2) and linearly with wind speed. Increasing the mid-canopy wind speed from 0 to 2 m s-1 changed the hydrometeor trajectory from -90° (vertical) to -26°. This resulted in a corresponding increase in the apparent canopy cover measured using voxel ray sampling from 0.2 to 0.7 (increase of 0.5) for the sparse forest plot and from 0.38 to 0.9 (increase of 0.52) for the dense forest plot.
The average canopy coverage of each forest plot measured from vertical (nadir) is also shown to be important determining how the apparent canopy structure, depending on trajectory angle changes across various trajectory angles. For example, the apparent canopy coverage of forest plot PWL_SW, which has the highest canopy coverage of 0.4 (-90°, nadir), shown in Figure 14 (b) increases at a faster rate with wind speed and reaches near complete apparent canopy coverage at lower wind speeds compared to sparser canopies. All forest plots approach completely closed canopies at horizontal trajectories (0°). However, this relationship does not consider the turbulent windflow patterns that exist within and above the canopy.
Discussion
The fraction of snowfall intercepted intercepted in this discontinuous subalpine forest is shown to be governed by forest structure, wind speed and initial canopy snow load. Occasionally the meteorology of a given event was observed to outweigh the differences in forest structure between different locations.
Satterlund & Haupt (1967) showed that interception efficiency increases as the canopy fills with snow bridging gaps in the canopy, while later declining due to branch bending and decreased canopy coverage and increased unloading. Hedstrom & Pomeroy (1998) and Storck et al. (2002) did not observe this initial increase. Observations from the subcanopy lysimeters here shown in Figure 8 (b) and Figure 9, despite displaying notable scatter, corroborate with the Satterlund & Haupt (1967) theory. Figure 8 (b) and Figure 9 indicate that initial canopy snow load exerted an important control on interception efficiency. A significant change of approximately 20% was detected between peak snow loads exceeding 20 mm (associated with low interception efficiency) and optimal snow loads of near 7 mm (associated with high interception efficiency). However, the decline in interception efficiency at high canopy snow loads in Figure 8 (b) and Figure 9 was rate much slower than has been observed by previous studies (Hedstrom & Pomeroy, 1998; Moeser et al., 2015; Satterlund & Haupt, 1967). At higher canopy snow loads, several studies suggest that interception efficiency declines (Hedstrom & Pomeroy, 1998; Moeser et al., 2015; Satterlund & Haupt, 1967), based on the premise of reduced canopy coverage due and catch efficiency to branch bending. However, the strong exponential decline in the interception efficiency observed with increasing event snowfall shown in Figure 1 may be a result of increased unloading rates as branches bend down. This may be attributed to the care taken in the observations presented here to minimize inclusion of unloading or ablation in the throughfall measurements. The longer duration between site visits in (Hedstrom & Pomeroy, 1998; Satterlund & Haupt, 1967) increase the likelihood of ablation in their throughfall measurements compared to the automated 15 minute measurements in this study. This potential inclusion of unloading within the interception parameterizations provided in (Hedstrom & Pomeroy, 1998; Moeser et al., 2015; Satterlund & Haupt, 1967) may lead to double counting of unloading when combined with an additional unloading parameterization as in (Hedstrom & Pomeroy, 1998). Still, if interception efficiency is expected to decline due to an associated decrease in canopy coverage from branch bending, as observed in Figure 9, it may be appropriate to have some decline in interception efficiency at higher canopy loads while also having a separate parameterization to increase unloading with canopy snow load.
Existing theories suggests either a positive (Storck et al., 2002) or negative (Hedstrom & Pomeroy, 1998) association of air temperature with the maximum canopy storage (Figure 2). The findings presented in this study deviate from the relationships outlined in Hedstrom & Pomeroy (1998) and Storck et al. (2002). Specifically, both Figure 8 (b) and Figure 9 showed an absence of a maximum canopy storage capacity. It is likely that a maximum canopy storage value does exist in the study site but was higher than the observed range of canopy storage values presented here. Additionally, in Figure 9, no relationship between air temperature and interception efficiency was observed. As maximum canopy storage value was not present in Figure 8 (a) and Figure 8 (b), a direct link cannot be made between air temperature and \(L_{max}\). However, the absence of a discernible relationship between interception efficiency and air temperature suggests that there might not be an influence on \(L_{max}\). This absence could also be explained by the concurrent presence of positive (adhesive and cohesive forces, Storck et al., 2002) and negative (branch bending, Schmidt & Pomeroy, 1990) relationships. These two opposing effects could balance each other out, resulting in no observable effect on interception efficiency in the observations presented in Figure 9. It is also plausible that the influence of wind on interception at this site masked any potential effect of air temperature.
Empirical evidence to support the theory of increased in canopy contact area with increased wind speed (Hedstrom & Pomeroy, 1998) has not been provided in the literature. Automated measurements of interception efficiency showed a slight increase in interception efficiency with increasing wind speed (below 1.5 m s-1). At wind speeds above 1.5 m s-1 interception efficiency was observed to decrease due to an associated increase in unloading. Measurements of interception efficiency collected by UAV-lidar surveys for a snowfall event with a 1.4 m s-1 average wind speed showed a large fraction of snow intercepted on the lee side of individual trees as a result of predominately non-vertical hydrometeor trajectory angles. The dominant hydrometeor trajectory angle was shown to be predicted using event wind speed and direction alone. The plot scale canopy coverage was shown to increase substantially as a result of increasing horizontal angle of hydrometers. The increase in plot scale canopy coverage observed in this study, of from 0.4 with zero wind to near complete coverage with midcanopy wind speeds of 2 m s-1 is much higher than the increase calculated using the method proposed in Hedstrom & Pomeroy (1998) to scale canopy coverage with wind speed. Using Equation 10 in Hedstrom & Pomeroy (1998), an increase of in canopy coverage of 0.03 is calculated using the same wind speed and forest structure and canopy height of 10 m and forested downwind distance of 100 m. The increase in interception efficiency observed away from individual trees in ?@fig-lidar-ip-site and the measured increase plot scale canopy coverage using voxel ray sampling in Figure 14 (b) show that Equation 10 in Hedstrom & Pomeroy (1998) may not be appropriate for sparse forests. Further investigation is required to test the relationships observed here in dense forest canopies. The large sensitivity of apparent forest cover to mid-canopy wind speed for this discontinuous subalpine forest suggests that this effect should be considered in snow interception models. While, increasing wind has been shown here to increase plot-scale interception, following the deposition of snow in the canopy wind is also known to increase rates of unloading (Bartlett & Verseghy, 2015; Roesch et al., 2001). Poor performance of the Roesch et al. (2001) and Bartlett & Verseghy (2015) algorithms has been shown in the literature (e.g. Lumbrazo et al., 2022). Since existing studies have not considered the concurrent effect of interception increasing as a result of increasing apparent forest structure due to increasing wind, validations of subsequent ablation parameterizations will likely benefit.
Conclusions
- Forest structure is the main factor in governing the fraction of intercepted snowfall at a particular site, with meteorological conditions contributing less to variability.
- Forest structure metrics calculated using high zenith angles better described the variability in interception efficiency in this windswept subalpine forest.
- Interception efficiency increased with increasing canopy coverage and LAI. However, the strength of this association was reduced at higher wind speeds.
- High wind speeds increased interception efficiency due to an associated increase in canopy contact area for sparsely forested locations, while later decreasing intercepted load due to increased snow unloading.
- No influence of air temperature on interception efficiency or maximum canopy storage was observed.
- Interception efficiency increased slightly as the canopy filled with snow and declined later at higher loads.
- Unloading increased with increasing canopy storage.
- The maximum canopy storage capacity for this study site is likely higher than was observed here.
References
Andreadis, K. M., Storck, P., & Lettenmaier, D. P. (2009). Modeling snow accumulation and ablation processes in forested environments. Water Resources Research, 45(5), 1–33. https://doi.org/10.1029/2008WR007042
Axelsson, P. (2000). DEM Generation from laser scanner data using adaptive TIN models. International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences - ISPRS Archives, 33, 110–117.
Bartlett, P. A., & Verseghy, D. L. (2015). Modified treatment of intercepted snow improves the simulated forest albedo in the Canadian Land Surface Scheme. Hydrological Processes, 29(14), 3208–3226. https://doi.org/10.1002/HYP.10431
Chen, J. M., Rich, P. M., Gower, S. T., Norman, J. M., & Plummer, S. (1997). Leaf area index of boreal forests: Theory, techniques, and measurements. Journal of Geophysical Research Atmospheres, 102(24), 29429–29443. https://doi.org/10.1029/97jd01107
Chianucci, F., & Macek, M. (2023). hemispheR: an R package for fisheye canopy image analysis. Agricultural and Forest Meteorology. https://doi.org/10.1016/j.agrformet.2023.109470
Cionco, R. M. (1965). A Mathematical Model for Air Flow in a Vegetative Canopy. Journal of Applied Meteorology (1962), 4(4), 517–522.
Cionco, R. M. (1972). A wind-profile index for canopy flow. Boundary-Layer Meteorology, 3(2), 255–263. https://doi.org/10.1007/BF02033923
Clark, M. P., Nijssen, B., Lundquist, J. D., Kavetski, D., Rupp, D. E., Woods, R. A., Freer, J. E., Gutmann, E. D., Wood, A. W., Brekke, L. D., Arnold, J. R., Gochis, D. J., & Rasmussen, R. M. (2015). A unified approach for process-based hydrologic modeling: 1. Modeling concept. Water Resources Research, 51(4), 2498–2514. https://doi.org/https://doi.org/10.1002/2015WR017198
Deems, J. S., Painter, T. H., & Finnegan, D. C. (2013). Lidar measurement of snow depth: A review. Journal of Glaciology, 59(215), 467–479. https://doi.org/10.3189/2013JoG12J154
Ellis, C. R., Pomeroy, J. W., Brown, T., & MacDonald, J. (2010). Simulation of snow accumulation and melt in needleleaf forest environments. Hydrology and Earth System Sciences, 14(6), 925–940. https://doi.org/10.5194/hess-14-925-2010
Ellis, C. R., Pomeroy, J. W., & Link, T. E. (2013). Modeling increases in snowmelt yield and desynchronization resulting from forest gap-thinning treatments in a northern mountain headwater basin. Water Resources Research, 49(2), 936–949. https://doi.org/10.1002/wrcr.20089
Essery, R., & Pomeroy, J. W. (2004). Vegetation and topographic control of wind-blown snow distributions in distributed and aggregated simulations for an arctic tundra basin. Journal of Hydrometeorology, 5(5), 735–744. https://doi.org/10.1175/1525-7541(2004)005<0735:VATCOW>2.0.CO;2
Essery, R., Pomeroy, J. W., Parviainen, J., & Storck, P. (2003). Sublimation of snow from coniferous forests in a climate model. Journal of Climate, 16(11), 1855–1864. https://doi.org/10.1175/1520-0442(2003)016<1855:SOSFCF>2.0.CO;2
Gelfan, A. N., Pomeroy, J. W., & Kuchment, L. S. (2004). Modeling forest cover influences on snow accumulation, sublimation, and melt. Journal of Hydrometeorology, 5(5), 785–803. https://doi.org/10.1175/1525-7541(2004)005<0785:MFCIOS>2.0.CO;2
Harder, P., & Pomeroy, J. W. (2013). Estimating precipitation phase using a psychrometric energy balance method. Hydrological Processes, 27(13), 1901–1914.
Harder, P., Pomeroy, J. W., & Helgason, W. D. (2020). Improving sub-canopy snow depth mapping with unmanned aerial vehicles: Lidar versus structure-from-motion techniques. Cryosphere, 14(6), 1919–1935. https://doi.org/10.5194/tc-14-1919-2020
Hedstrom, N. R., & Pomeroy, J. W. (1998). Measurements and modelling of snow interception in the boreal forest. Hydrological Processes, 12(10-11), 1611–1625. https://doi.org/10.1002/(SICI)1099-1085(199808/09)12:10/11<1611::AID-HYP684>3.0.CO;2-4
Helbig, N., Moeser, D., Teich, M., Vincent, L., Lejeune, Y., Sicart, J. E. J.-E. J. E., & Monnet, J. M. J.-M. (2020). Snow processes in mountain forests: interception modeling for coarse-scale applications. Hydrology and Earth System Sciences, 24(5), 2545–2560. https://doi.org/10.5194/hess-24-2545-2020
Hijmans, R. J. (2024). terra: Spatial Data Analysis. https://cran.r-project.org/package=terra
Katsushima, T., Kato, A., Aiura, H., Nanko, K., Suzuki, S., Takeuchi, Y., & Murakami, S. (2023). Modelling of snow interception on a Japanese cedar canopy based on weighing tree experiment in a warm winter region. Hydrological Processes, 37(6), 1–16. https://doi.org/10.1002/hyp.14922
Krinner, G., Derksen, C., Essery, R., Flanner, M., Hagemann, S., Clark, M. P., Hall, A., Rott, H., Brutel-Vuilmet, C., Kim, H., Ménard, C. B., Mudryk, L., Thackeray, C., Wang, L., Arduini, G., Balsamo, G., Bartlett, P., Boike, J., Boone, A., … Zhu, D. (2018). ESM-SnowMIP: Assessing snow models and quantifying snow-related climate feedbacks. Geoscientific Model Development, 11(12), 5027–5049. https://doi.org/10.5194/gmd-11-5027-2018
Langs, L. E., Petrone, R. M., & Pomeroy, J. W. (2020). A \(\delta\)18O and \(\delta\)2H stable water isotope analysis of subalpine forest water sources under seasonal and hydrological stress in the Canadian Rocky Mountains. Hydrological Processes, 34(26), 5642–5658. https://doi.org/10.1002/hyp.13986
LAStools. (2024). Efficient LiDAR Processing Software (version 240220, academic). http://rapidlasso.com/LAStools
Lumbrazo, C., Bennett, A., Currier, W. R., Nijssen, B., & Lundquist, J. (2022). Evaluating Multiple Canopy-Snow Unloading Parameterizations in SUMMA With Time-Lapse Photography Characterized by Citizen Scientists. Water Resources Research, 58(6), 1–22. https://doi.org/10.1029/2021WR030852
Lundquist, J. D., Dickerson-Lange, S., Gutmann, E., Jonas, T., Lumbrazo, C., & Reynolds, D. (2021). Snow interception modelling: Isolated observations have led to many land surface models lacking appropriate temperature sensitivities. Hydrological Processes, 35(7), 1–20. https://doi.org/10.1002/hyp.14274
Macdonald, J., & Pomeroy, J. W. (2007). Gauge Undercatch of Two Common Snowfall Gauges in a Prairie Environment. Proceedings of the 64th Eastern Snow Conference, St. John‘s, Canada., 1974, 119–126.
Moeser, D., Stähli, M., & Jonas, T. (2015). Improved snow interception modeling using canopy parameters derived from airborne LiDAR data. Water Resources Research, 51(7), 5041–5059.
Natural Resources Canada. (2024). Precise point positioning. https://webapp.csrs-scrs.nrcan-rncan.gc.ca/geod/tools-outils/ppp.php
Parviainen, J., & Pomeroy, J. W. (2000). Multiple-scale modelling of forest snow sublimation: Initial findings. Hydrological Processes, 14(15), 2669–2681. https://doi.org/10.1002/1099-1085(20001030)14:15<2669::AID-HYP85>3.0.CO;2-Q
Pomeroy, J. W., Parviainen, J., Hedstrom, N., & Gray, D. M. (1998). Coupled modelling of forest snow interception and sublimation. Hydrological Processes, 12(15), 2317–2337. https://doi.org/10.1002/(SICI)1099-1085(199812)12:15<2317::AID-HYP799>3.0.CO;2-X
Rasouli, K., Pomeroy, J. W., & Whitfield, P. H. (2019). Are the effects of vegetation and soil changes as important as climate change impacts on hydrological processes? Hydrology and Earth System Sciences, 23(12), 4933–4954. https://doi.org/10.5194/hess-23-4933-2019
Roesch, A., Wild, M., Gilgen, H., & Ohmura, A. (2001). A new snow cover fraction parameterization for the ECHAM4 GCM. Climate Dynamics, 17(12), 933–946. https://doi.org/10.1007/s003820100153
Roth, T. R., & Nolin, A. W. (2019). Characterizing Maritime Snow Canopy Interception in Forested Mountains. Water Resources Research, 55(6), 4564–4581. https://doi.org/10.1029/2018WR024089
Rutter, N., Essery, R., Pomeroy, J. W., Altimir, N., Andreadis, K. M., Baker, I., Barr, A., Bartlett, P., Boone, A., Deng, H., Douville, H., Dutra, E., Elder, K., Ellis, C. R., Feng, X., Gelfan, A., Goodbody, A., Gusev, Y., Gustafsson, D., … Yamazaki, T. (2009). Evaluation of forest snow processes models (SnowMIP2). Journal of Geophysical Research: Atmospheres, 114(D6), 10–18. https://doi.org/10.1029/2008JD011063
Satterlund, D. R., & Haupt, H. F. (1967). Snow catch by Conifer Crowns. Water Resources Research, 3(4), 1035–1039.
Schmidt, R. A., & Pomeroy, J. W. (1990). Bending of a conifer branch at subfreezing temperatures: implications for snow interception. Canadian Journal of Forest Research, 20(8), 1251–1253. https://doi.org/10.1139/x90-165
Staines, J., & Pomeroy, J. W. (2023). Influence of forest canopy structure and wind flow on patterns of sub-canopy snow accumulation in montane needleleaf forests. Hydrological Processes, 37(10), 1–19. https://doi.org/10.1002/hyp.15005
Storck, P., Lettenmaier, D. P., & Bolton, S. M. (2002). Measurement of snow interception and canopy effects on snow accumulation and melt in a mountainous maritime climate, Oregon, United States. Water Resources Research, 38(11), 1–16. https://doi.org/10.1029/2002wr001281
Zhu, J., Matsuzaki, T., & Gonda, Y. (2001). Wind profiles in a coastal forest of Japanese black pine (Pinus thunbergii Parl.) with different thinning intensities. Journal of Forest Research, 6(4), 287–296. https://doi.org/10.1007/BF02762470